/*
// Description
// ----------------------------
This program performs a spatial RD analysis on merged Experian 
EV sales data for Census block groups on the borders of IOU
and municipal utilities in California from 2014-2017. */


/*** START CODE ***/


// Preamble
// ----------------------------
capture log close
clear all
*version 14.2
set more off
set matsize 800

set scheme s1color

use "${Data_Clean}/spatial_RD_regdata.dta", clear

local delta_demog "delta_pop_density delta_mud_hh_share "
local delta_dist "cons i.groupid#c.IOU_dist i.groupid#c.muni_dist"

forvalues x = 2/5 {
	gen IOU_dist`x' = IOU_dist^`x'
	gen muni_dist`x' = muni_dist^`x'
}

local other_cov "delta_fuel_economy_mean delta_hybrid_share delta_luxury_share delta_income delta_population" 
local iou_fe "PGE_dum SCE_dum SDGE_dum"

local dist_km = 8
local dist_m = `dist_km'*1000
local dist_condition "distance_btw_cbg0 < `dist_m'"

local gas_price_var "delta_obswgtprice_lt3_"
local elec_price_var "delta_p_e_active"
local price_var "`elec_price_var' `gas_price_var'"

local icempg = 30
local BEVmpkwh = 4

gen delta_p_e_active_permile = `elec_price_var' / `BEVmpkwh'
gen delta_P_gas_permile = `gas_price_var' / `icempg'

label var delta_p_e_active_permile "$\Delta$ Marg. Price (cents/mile)"
label var delta_P_gas_permile "$\Delta$ Gas Price (cents/mile)"

local price_var_permile "delta_p_e_active_permile delta_P_gas_permile"

estimates clear

/* Main RD Results Table*/

reghdfe delta_share_BEV `price_var', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj1

reghdfe delta_share_BEV `price_var' `delta_demog' , absorb(`delta_dist')  vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj2

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj3

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' `iou_fe', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estadd local IOU_FE "X" , replace
estimates store ref_and_adj4

/* Main RD Results Table - within 8k of boundary*/

reghdfe delta_share_BEV `price_var' if `dist_condition', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj1_close

reghdfe delta_share_BEV `price_var' `delta_demog' if `dist_condition' , absorb(`delta_dist')  vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj2_close

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if `dist_condition', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj3_close

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' `iou_fe' if `dist_condition', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estadd local IOU_FE "X" , replace
estimates store ref_and_adj4_close


/* Main RD Results Table - Per Mile Gas / Electricity Costs*/

reghdfe delta_share_BEV `price_var_permile', absorb(`delta_dist') vce(cluster cbg0 cbg1)
test delta_p_e_active_permile=delta_P_gas_permile
estadd local testp = round(r(p),.001) , replace
nlcom -_b[delta_p_e_active_permile]/_b[delta_P_gas_permile]
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj1_permile

reghdfe delta_share_BEV `price_var_permile' `delta_demog' , absorb(`delta_dist')  vce(cluster cbg0 cbg1)
test delta_p_e_active_permile=delta_P_gas_permile
estadd local testp = round(r(p),.001) , replace
nlcom -_b[delta_p_e_active_permile]/_b[delta_P_gas_permile]
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj2_permile

reghdfe delta_share_BEV `price_var_permile' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
test delta_p_e_active_permile=delta_P_gas_permile
estadd local testp = round(r(p),.001) , replace
nlcom -_b[delta_p_e_active_permile]/_b[delta_P_gas_permile]
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj3_permile

reghdfe delta_share_BEV `price_var_permile' `delta_demog' `other_cov' `iou_fe', absorb(`delta_dist') vce(cluster cbg0 cbg1)
test delta_p_e_active_permile=delta_P_gas_permile
estadd local testp = round(r(p),.001) , replace
nlcom -_b[delta_p_e_active_permile]/_b[delta_P_gas_permile]
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estadd local IOU_FE "X" , replace
estimates store ref_and_adj4_permile


/* Robustness check table*/

/* Outside PGE */
reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if PGE_dum==0, absorb(`delta_dist')  vce(cluster cbg0 cbg1)
estadd local dist_cont "Exc. PGE" , replace
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj3a

*Graphs of differential Population Density and MUD by MUNI / IOU
preserve

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
gen insample=e(sample)
keep if insample==1 

gen IOU_name = ""
replace IOU_name = utility_name0 if ref_IOU==1
replace IOU_name = utility_name1 if ref_IOU==0
gen Muni_name = ""
replace Muni_name = utility_name1 if ref_IOU==1
replace Muni_name = utility_name0 if ref_IOU==0
tab IOU_name
tab Muni_name

*Define differentials where differentials are Muni Demographic - IOU demographic
gen mud_hh_diff = .
replace mud_hh_diff = delta_mud_hh_share if ref_IOU==0
replace mud_hh_diff = -delta_mud_hh_share if ref_IOU==1

gen pop_dens_diff = .
replace pop_dens_diff = delta_pop_density if ref_IOU==0
replace pop_dens_diff = -delta_pop_density if ref_IOU==1

gen totalpop = population0+population1

collapse (mean) pop_dens_diff mud_hh_diff, by(Muni_name)
sum pop_dens, det
scalar PDcut = r(p75)
sum mud_hh_diff, det
scalar mudcut = r(p75)
macro list
scalar list
gen density_subset = (mud_hh_diff<mudcut) & (pop_dens_diff<PDcut)

/*Figure A5*/
graph twoway (scatter pop_dens_diff mud_hh_diff if density_subset==0, mlab(Muni_name) mlabpos(1) mlabsize(vsmall) msize(small) msymb(Dh) color(gs3)) ///
(scatter pop_dens_diff mud_hh_diff if density_subset==1, mlab(Muni_name) mlabpos(11) mlabsize(vsmall) msize(small) msymb(Dh) color(olive)) ///
(scatteri `=PDcut' -.45 `=PDcut' `=mudcut' -2.5 `=mudcut', lp(dash) c(l) lc(gs8%50) m(i)), ///
ytitle("Difference in Population Density") xtitle("Difference in Share of Multi-unit Dwellings") ///
xscale(range(-.5 .5)) legend(off)
graph export "$ResultsOut/Figures/RD_subsample.png", replace

tempfile bordersubset
keep Muni_name density_subset
save `bordersubset', replace

restore

gen Muni_name = ""
replace Muni_name = utility_name1 if ref_IOU==1
replace Muni_name = utility_name0 if ref_IOU==0

merge n:1 Muni_name using `bordersubset'
tab _merge
drop _merge

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if density_subset==1, absorb(`delta_dist')  vce(cluster cbg0 cbg1)
estadd local dist_cont "Exc. PGE" , replace
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj3b  


/*Generate regression tables*/

/*Table 2*/	
local reportvars "`price_var' `delta_demog' `other_cov'"

esttab ref_and_adj1 ref_and_adj2 ref_and_adj3 ref_and_adj4 ref_and_adj1_close ref_and_adj2_close ref_and_adj3_close ref_and_adj4_close ///
	using "${ResultsOut}/Tables/RD_results.tex", ///
	keep(`reportvars') ///
	replace label se star(* 0.10 ** 0.05 *** 0.01) nogaps nonotes noconstant booktabs ///
	b(a2) se(a2) stats(IOU_FE impliedgamma gammase N r2,  labels("IOU FE"  "\\ Implied $\gamma$"   " " "\\ Observations" "R-Squared")) ///
    mgroups("\underline{Full Sample}" "\underline{CBG dist $<$  `dist_km'km}", pattern(1 0 0 0 1 0 0 0 ) span prefix(\multicolumn{@span}{c}{) suffix(}))  ///
	nomtitles 
		
/*Table 3*/	
local reportvars "`price_var' `delta_demog' `other_cov'"

esttab ref_and_adj3 ref_and_adj3a ref_and_adj3b ///
	using "${ResultsOut}/Tables/RD_results_both_robustness.tex", ///
	keep(`reportvars') ///
	replace label se star(* 0.10 ** 0.05 *** 0.01) nogaps nonotes noconstant booktabs ///
	b(a2) se(a2) stats(impliedgamma gammase N r2,  labels("Implied $\gamma$"   " " "\\ Observations" "R-Squared")) ///
	mtitles("Full Sample" "Excl. PGE CBGs" "Exc. Urban/Rural Boundaries")
	
	 
	 
/*Main RD Results Table PHEVs*/

reghdfe delta_share_PHEV `price_var', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj1_PHEV

reghdfe delta_share_PHEV `price_var' `delta_demog' , absorb(`delta_dist')  vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj2_PHEV

reghdfe delta_share_PHEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj3_PHEV

reghdfe delta_share_PHEV `price_var' `delta_demog' `other_cov' `iou_fe', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estadd local IOU_FE "X" , replace
estimates store ref_and_adj4_PHEV


/* Main RD Results Table - within 8k of boundary*/

reghdfe delta_share_PHEV `price_var' if `dist_condition', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj1_close_PHEV

reghdfe delta_share_PHEV `price_var' `delta_demog' if `dist_condition' , absorb(`delta_dist')  vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj2_close_PHEV

reghdfe delta_share_PHEV `price_var' `delta_demog' `other_cov' if `dist_condition', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store ref_and_adj3_close_PHEV

reghdfe delta_share_PHEV `price_var' `delta_demog' `other_cov' `iou_fe' if `dist_condition', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estadd local IOU_FE "X" , replace
estimates store ref_and_adj4_close_PHEV
	

local reportvars "`price_var' `delta_demog' `other_cov'"

esttab ref_and_adj1_PHEV ref_and_adj2_PHEV ref_and_adj3_PHEV ref_and_adj4_PHEV ref_and_adj1_close_PHEV ref_and_adj2_close_PHEV ref_and_adj3_close_PHEV ref_and_adj4_close_PHEV ///
	using "${ResultsOut}/Tables/RD_results_PHEV.tex", ///
	keep(`reportvars') ///
	replace label se star(* 0.10 ** 0.05 *** 0.01) nogaps nonotes noconstant booktabs ///
	b(a2) se(a2) stats(IOU_FE impliedgamma gammase N r2,  labels("IOU FE"  "\\ Implied $\gamma$"   " " "\\ Observations" "R-Squared")) ///
    mgroups("\underline{Full Sample}" "\underline{CBG dist $<$  `dist_km'km}", pattern(1 0 0 0 1 0 0 0 ) span prefix(\multicolumn{@span}{c}{) suffix(}))  ///
	nomtitles 
	
		 
/*Figure 5, Appendix Figure A3*/
/*Graphs of different bandwidths */

preserve
clear
set obs 20
gen dist_bw = _n
gen gamma = .
gen gamma_se = .
gen elec_coef = .
gen elec_coef_se = .
gen gas_coef = .
gen gas_coef_se = .
tempfile bwgraph
save `bwgraph', replace
restore

forvalues x = 1/20 {
	preserve
	local dist_bw = `x'*1000
local dist_bw_condition "distance_btw_cbg0 < `dist_bw'"
reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if `dist_bw_condition', absorb(`delta_dist') vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')
local impliedgamma=r(b)[1,1]
local gammase=sqrt(r(V)[1,1])
macro list

use `bwgraph', clear
replace gamma = `impliedgamma' if dist_bw==`x' 
replace gamma_se = `gammase' if dist_bw==`x'
replace elec_coef = _b[delta_p_e_active] if dist_bw==`x'
replace elec_coef_se = _se[delta_p_e_active] if dist_bw==`x'
replace gas_coef = _b[`gas_price_var'] if dist_bw==`x'
replace gas_coef_se = _se[`gas_price_var'] if dist_bw==`x'
save `bwgraph', replace
restore

}	

preserve
use `bwgraph',clear
foreach x in  elec_coef gas_coef gamma {
	gen `x'_topCI = `x' + 1.96*`x'_se
	gen `x'_bottomCI = `x' - 1.96*`x'_se
}

/*Figure 5*/
graph twoway (rarea gamma_topCI gamma_bottomCI dist_bw, col(gs8%40) lw(none)) /// 
(line gamma dist_bw, lc(black)), ///
xtitle("Distance Bandwidth (km)") ytitle("Implied Gamma") ///
legend(off) yline(0 1, lp(dash) lc(gs2%90)) yscale(range(-.4 1)) ylabel(-.4(.2)1.2)
graph export "${ResultsOut}/Figures/gamma_RDbandwidth.png", replace

/*Figure A3*/
graph twoway (rarea elec_coef_topCI elec_coef_bottomCI dist_bw, col(gs8%40) lw(none)) /// 
(line elec_coef dist_bw, lc(black)), ///
xtitle("Distance Bandwidth (km)") ytitle("Electricity Price Coefficient") ///
legend(off) yline(0, lp(dash))
graph export "${ResultsOut}/Figures/elec_coef_RDbandwidth.png", replace

/*Figure A3*/
graph twoway (rarea gas_coef_topCI gas_coef_bottomCI dist_bw, col(gs8%40) lw(none)) ///
(line gas_coef dist_bw, lc(black)), ///
xtitle("Distance Bandwidth (km)") ytitle("Gasoline Price Coefficient") ///
legend(off) yline(0, lp(dash))
graph export "${ResultsOut}/Figures/gas_coef_RDbandwidth.png", replace


restore
	
/*Figures 6, Appendix Figure A4, Appendix Table A1*/	
/*Robustness checks for alternative polynomials */

estimates clear

	local delta_dist1 "cons i.groupid#c.IOU_dist i.groupid#c.muni_dist"
	local delta_dist2 "`delta_dist1' i.groupid#c.IOU_dist2 i.groupid#c.muni_dist2"
	local delta_dist3 "`delta_dist2' i.groupid#c.IOU_dist3 i.groupid#c.muni_dist3"
	local delta_dist4 "`delta_dist3' i.groupid#c.IOU_dist4 i.groupid#c.muni_dist4"
	local delta_dist5 "`delta_dist4' i.groupid#c.IOU_dist5 i.groupid#c.muni_dist5"	
	
	macro list
		
/*Baseline Main Specification*/

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist1') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly1coef
nlcom (nl_1:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_1:_b[delta_p_e_active]) (pg_1:_b[`gas_price_var']), post
gen insample=e(sample)
estimates store RD_robust_poly1

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if `dist_condition', absorb(`delta_dist1') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly_near1coef
nlcom (nl_1:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_1:_b[delta_p_e_active]) (pg_1:_b[`gas_price_var']), post
gen insample_dist=e(sample)
estimates store RD_robust_poly_near1


reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if insample==1, absorb(groupid) vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly0coef
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_0:_b[delta_p_e_active]) (pg_0:_b[`gas_price_var']), post
estimates store RD_robust_poly0

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if `dist_condition' & insample_dist==1, absorb(groupid) vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly_near0coef
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_0:_b[delta_p_e_active]) (pg_0:_b[`gas_price_var']), post
estimates store RD_robust_poly_near0


reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist2') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly2coef
nlcom (nl_2:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_2:_b[delta_p_e_active]) (pg_2:_b[`gas_price_var']), post
estimates store RD_robust_poly2

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if `dist_condition', absorb(`delta_dist2') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly_near2coef
nlcom (nl_2:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_2:_b[delta_p_e_active]) (pg_2:_b[`gas_price_var']), post
estimates store RD_robust_poly_near2


reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist3') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly3coef
nlcom (nl_3:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_3:_b[delta_p_e_active]) (pg_3:_b[`gas_price_var']), post
estimates store RD_robust_poly3

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if `dist_condition', absorb(`delta_dist3') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly_near3coef
nlcom (nl_3:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_3:_b[delta_p_e_active]) (pg_3:_b[`gas_price_var']), post
estimates store RD_robust_poly_near3


reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist4') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly4coef
nlcom (nl_4:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_4:_b[delta_p_e_active]) (pg_4:_b[`gas_price_var']), post
estimates store RD_robust_poly4

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov' if `dist_condition', absorb(`delta_dist4') vce(cluster cbg0 cbg1)
nlcom (nl_0:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) 
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store RD_robust_poly_near4coef
nlcom (nl_4:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')) (pe_4:_b[delta_p_e_active]) (pg_4:_b[`gas_price_var']), post
estimates store RD_robust_poly_near4

/*Figure 6*/
coefplot (RD_robust_poly0 RD_robust_poly1 RD_robust_poly2 RD_robust_poly3 RD_robust_poly4) /// 
(RD_robust_poly_near0 RD_robust_poly_near1  RD_robust_poly_near2  RD_robust_poly_near3  RD_robust_poly_near4), /// 
keep(nl_0 nl_1 nl_2 nl_3 nl_4) vertical yscale(range(-0.2 1.2)) ylabel(-0.2 (0.2) 1.2) yline(0 1, lp(dash) lc(gs2))  ///
coeflabels(nl_0 = "Distance Exc." nl_1 = "Linear" nl_2 = "2nd Order" nl_3 = "3rd Order" nl_4 = "4th Order") ///
legend(order(2 "Full Sample" 4 "Within 8km")) ytitle("Implied Gamma") 

graph export "${ResultsOut}/Figures/gamma_RDpolynomial.png", replace


/*Figure A4*/
coefplot (RD_robust_poly0 RD_robust_poly1 RD_robust_poly2 RD_robust_poly3 RD_robust_poly4) /// 
(RD_robust_poly_near0 RD_robust_poly_near1  RD_robust_poly_near2  RD_robust_poly_near3  RD_robust_poly_near4), /// 
keep(pe_0 pe_1 pe_2 pe_3 pe_4) vertical yscale(range(-0.6 0.6)) ylabel(-0.6 (0.2) 0.6) yline(0 1, lp(dash) lc(gs2))  ///
coeflabels(pe_0 = "Distance Exc." pe_1 = "Linear" pe_2 = "2nd Order" pe_3 = "3rd Order" pe_4 = "4th Order") ///
legend(order(2 "Full Sample" 4 "Within 8km")) ytitle("Electricity Price Coefficient") 

graph export "${ResultsOut}/Figures/elec_coef_RDpolynomial.png", replace

coefplot (RD_robust_poly0 RD_robust_poly1 RD_robust_poly2 RD_robust_poly3 RD_robust_poly4) /// 
(RD_robust_poly_near0 RD_robust_poly_near1  RD_robust_poly_near2  RD_robust_poly_near3  RD_robust_poly_near4), /// 
keep(pg_0 pg_1 pg_2 pg_3 pg_4) vertical yscale(range(-0.6 0.6)) ylabel(-0.6 (0.2) 0.6) yline(0 1, lp(dash) lc(gs2))   ///
coeflabels(pg_0 = "Distance Exc." pg_1 = "Linear" pg_2 = "2nd Order" pg_3 = "3rd Order" pg_4 = "4th Order") ///
legend(order(2 "Full Sample" 4 "Within 8km")) ytitle("Gasoline Price Coefficient") 

graph export "${ResultsOut}/Figures/gas_coef_RDpolynomial.png", replace

/*Table A1*/	
local reportvars "`price_var' `delta_demog' `other_cov'"

esttab RD_robust_poly0coef RD_robust_poly1coef RD_robust_poly2coef RD_robust_poly3coef RD_robust_poly4coef RD_robust_poly_near0coef RD_robust_poly_near1coef RD_robust_poly_near2coef RD_robust_poly_near3coef RD_robust_poly_near4coef  ///
	using "${ResultsOut}/Tables/RD_robustness_poly.tex", ///
	keep(`reportvars' ) ///
	rename(`renamevars') ///
	replace label se star(* 0.10 ** 0.05 *** 0.01) nogaps nonotes noconstant booktabs ///
		b(a2) se(a2) stats(impliedgamma gammase N r2,  labels("Implied $\gamma$"   " " "\\ Observations" "R-Squared")) ///
	mtitles("Dist Exc." "Linear" "2nd Degree" "3rd Degree" "4th Degree" "Dist Exc." "Linear" "2nd Degree" "3rd Degree" "4th Degree")  ///
    mgroups("\underline{Full RD Sample}" "\underline{CBG dist $<$  `dist_km'km}", pattern(1 0 0 0 0 1 0 0 0 0) span prefix(\multicolumn{@span}{c}{) suffix(})) 

	
/*Figure 7*/
/*Robustness checks using alternative assumptions*/

set scheme s1color
use "${Data_Clean}/spatial_RD_regdata.dta", clear

local delta_demog "delta_pop_density delta_mud_hh_share "
*local delta_dist "cons"
local delta_dist "cons i.groupid#c.IOU_dist i.groupid#c.muni_dist"

local other_cov "delta_fuel_economy_mean delta_hybrid_share delta_luxury_share delta_income delta_population" 
local iou_fe "PGE_dum SCE_dum SDGE_dum"

local dist_km = 8
local dist_m = `dist_km'*1000
local dist_condition "distance_btw_cbg0 < `dist_m'"

local gas_price_var "delta_obswgtprice_lt3_"
local elec_price_var "delta_p_e_active"
local price_var "`elec_price_var' `gas_price_var'"

local icempg = 30
local BEVmpkwh = 4

gen delta_p_e_active_permile = `elec_price_var' / `BEVmpkwh'
gen delta_P_gas_permile = `gas_price_var' / `icempg'

label var delta_p_e_active_permile "$\Delta$ Marg. Price (cents/mile)"
label var delta_P_gas_permile "$\Delta$ Gas Price (cents/mile)"

local price_var_permile "delta_p_e_active_permile delta_P_gas_permile"

estimates clear

/*Baseline Specification*/

reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_base
nlcom (gamma_base:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')), post
estimates store RD_robust_base

/*Alternative Replacement Vehicle*/
local alticempg = 37.1
reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_rep_veh_MR
nlcom (gamma_rep_veh_MR:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`alticempg')), post
estimates store RD_robust_rep_veh_MR

/*VMT Specification*/
/*VMT X% Lower with EV - Davis (2019)*/
/*EVMT BEVs 6300/10200, EVMT ALL EVs 7200 / 9800 =73.5%*/
local evmt = .618 /*Lower with EV - Davis (2019)*/
reghdfe delta_share_BEV delta_p_e_active `gas_price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_evmt
nlcom (gamma_evmt:-(_b[delta_p_e_active]*`BEVmpkwh'/`evmt')/(_b[`gas_price_var']*`icempg')), post
estimates store RD_robust_evmt

/*Electricity Price Robustness Checks*/
/*Electricity Prices - 20% away from home charging, free*/
local percenthome = .8 /*Citation from Hardiman et al. (2018)*/
gen delta_p_e_alt1 = `percenthome'*delta_p_e_active
reghdfe delta_share_BEV delta_p_e_alt1 `gas_price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_p_e_alt1
nlcom (gamma_p_e_alt1:-(_b[delta_p_e_alt1]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')), post
estimates store RD_robust_p_e_alt1

/*Electricity Prices - Average Price rather than marginal*/
/*Average consumption for EV owning household*/
gen delta_p_e_alt2 = delta_KWH900AP
reghdfe delta_share_BEV delta_p_e_alt2 `gas_price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_p_e_alt2
nlcom (gamma_p_e_alt2:-(_b[delta_p_e_alt2]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')), post
estimates store RD_robust_p_e_alt2
gen insample = e(sample)

/*Electricity Prices - Alternative Marginal Prices for IOUs, down a tier*/
gen delta_p_e_alt3 = 100*(AvgRateTier20 - p_e_active1) if ref_IOU==1
replace delta_p_e_alt3 = 100*(AvgRateTier00 -p_e_active1) if ref_IOU==1 & year>2016
replace delta_p_e_alt3 = 100*(p_e_active0 - AvgRateTier21) if ref_IOU==0
replace delta_p_e_alt3 = 100*(p_e_active0 - AvgRateTier01) if ref_IOU==0 & year>2016
reghdfe delta_share_BEV delta_p_e_alt3 `gas_price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_p_e_alt3
nlcom (gamma_p_e_alt3:-(_b[delta_p_e_alt3]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')), post
estimates store RD_robust_p_e_alt3

/*Electricity Prices - Alternative Marginal Prices for IOUs, use penalty tier*/
gen delta_p_e_alt4 = delta_p_e_high
reghdfe delta_share_BEV delta_p_e_alt4 `gas_price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_p_e_alt4
nlcom (gamma_p_e_alt4:-(_b[delta_p_e_alt4]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg')), post
estimates store RD_robust_p_e_alt4

/*Gasoline Prices - Robustness checks*/
/*50% shopping in lower price location*/
local shopping = .5
gen delta_obswgtprice_alt1 = `shopping'*delta_obswgtprice_lt3_
reghdfe delta_share_BEV delta_p_e_active delta_obswgtprice_alt1 `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_obswgtprice_alt1
nlcom (gamma_obswgtprice_alt1:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[delta_obswgtprice_alt1]*`icempg')), post
estimates store RD_obswgtprice_alt1

/*Asymmetry in Ownership Lengths*/
/*Scrappage 25% higher for EVs or lower for EVs*/

/*25% higher scrappage rate*/
local EVscrappagerate = 1.25
reghdfe delta_share_BEV `price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_scrappage1
nlcom (gamma_scrappage1:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg'*(`EVscrappagerate'))), post
estimates store RD_scrappage1

/*25% lower scrappage rate*/
local EVscrappagerate = 1/1.25 
reghdfe delta_share_BEV  `price_var' `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_scrappage2
nlcom (gamma_scrappage2:-(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`gas_price_var']*`icempg'*(`EVscrappagerate'))), post
estimates store RD_scrappage2

/*MULTIPLE*/
/*Electricity Prices - 20% away from home charging, free*/
/*VMT X% Lower with EV - Davis (2019)*/
/*Alternative Replacement Vehicle*/
reghdfe delta_share_BEV delta_p_e_alt1 delta_obswgtprice_lt3_ `delta_demog' `other_cov', absorb(`delta_dist') vce(cluster cbg0 cbg1)
estimates store RD_robust_p_e_multi
nlcom (gamma_p_e_multi:-(_b[delta_p_e_alt1]*`BEVmpkwh'/`evmt')/(_b[`gas_price_var']*`alticempg')), post
estimates store RD_robust_p_e_multi

/*Figure 7*/
coefplot (RD_robust_base RD_robust_rep_veh_MR RD_robust_evmt RD_robust_p_e_alt1 RD_robust_p_e_alt2 RD_robust_p_e_alt3 RD_robust_p_e_alt4 RD_obswgtprice_alt1 RD_scrappage1 RD_scrappage2 RD_robust_p_e_multi) , /// 
keep(gamma_base gamma_rep_veh_MR gamma_evmt gamma_p_e_alt1 gamma_p_e_alt2 gamma_p_e_alt3 gamma_p_e_alt4 gamma_obswgtprice_alt1 gamma_scrappage1 gamma_scrappage2 gamma_p_e_multi) /// 
coeflabels(gamma_base = "[1] Base Model" gamma_rep_veh_MR = "[2] Replacement MPG of 37.1" /// 
gamma_evmt = "[3] 38% Lower eVMT" /// 
gamma_p_e_alt1 = "[4] 20% Free Public Charging" /// 
gamma_p_e_alt2 = "[5] Average Residential Prices" ///
gamma_p_e_alt3 = "[6] IOU Lower Tier Prices" ///
gamma_p_e_alt4 = "[7]  IOU Penalty Tier Prices" ///
gamma_obswgtprice_alt1 = "[8] 50% Gasoline Shopping" /// 
gamma_scrappage1 = "[9] 25% Higher EV Scrappage" /// 
gamma_scrappage2 = "[10] 25% Lower EV Scrappage" ///
gamma_p_e_multi = "[11] Scenario [2] + [3] + [4]]" ///
) ///
headings(gamma_base = " "  gamma_rep_veh_MR = " " gamma_evmt = " " gamma_p_e_alt1 = " " gamma_obswgtprice_alt1 = " " gamma_scrappage1 = " " gamma_p_e_multi = " ", offset(0.25) nogap)  ///
horizontal xscale(range(-0.2 1.2)) xlabel(-0.2 (0.2) 1.2) xline(1, lp(dash) lc(gs2)) xline(0.157, lp(dash) lc(gs2%40)) ///
legend(off) xtitle("Implied Gamma") yscale(range(0 20)) scale(0.8)

graph export "${ResultsOut}/Figures/gamma_alt_assumptions.png", replace


	
/*Table A3*/	
/*Robustness checks for gas alternative prices*/
local ranges "lt1 lt3 lt5 lt10"

foreach x in `ranges' {
	rename delta_obswgtprice_`x'_ delta_P_ow_`x'
	local altgaspricelist "`altgaspricelist' delta_P_ow_`x'"
}

foreach x in `ranges' {
	rename delta_invdistwgtprice_`x'_ delta_P_idw_`x'
	local altgaspricelist "`altgaspricelist' delta_P_idw_`x'"
}

local fullgasprice_vars "delta_P_gas `altgaspricelist'"

local resultscovlist ""

foreach x of varlist `fullgasprice_vars' {
	
reghdfe delta_share_BEV delta_p_e_active `x' `delta_demog' `other_cov', absorb(`delta_dist')  vce(cluster cbg0 cbg1)
nlcom -(_b[delta_p_e_active]*`BEVmpkwh')/(_b[`x']*`icempg')
estadd local impliedgamma=round(r(b)[1,1],.001) , replace
local gammase=round(sqrt(r(V)[1,1]),.001)
estadd local gammase="("+"`gammase'"+")", replace
estimates store cov_altgas_`x'
local resultscovlist "`resultscovlist' cov_altgas_`x'"
}

local renamevars ""
foreach x of varlist `altgaspricelist' {
local renamevars " `renamevars' `x' delta_P_gas"
}

/*Table A3*/	
esttab `resultscovlist' ///
	using "${ResultsOut}/Tables/RD_robustness_gasprice_vars_cov.tex", ///
	keep(delta_p_e_active delta_P_gas `delta_demog'  `other_cov' ) ///
	rename(`renamevars') ///
	replace label se star(* 0.10 ** 0.05 *** 0.01) nogaps nonotes noconstant booktabs ///
		b(a2) se(a2) stats(impliedgamma gammase N r2,  labels("Implied $\gamma$"   " " "\\ Observations" "R-Squared")) ///
	mtitles("Avg Zip Price" "1mi Radius" "3mi Radius" "5mi Radius" "10mi Radius" "1mi Radius" "3mi Radius" "5mi Radius" "10mi Radius" )  ///
    mgroups("" "\underline{Observation Weighted Average Prices}" "\underline{Inverse Distance Weighted Average Prices}", pattern(1 1 0 0 0 1 0 0 0) span prefix(\multicolumn{@span}{c}{) suffix(})) 

	